#tfile="/home1/99/jc152199/MicroclimateStatisticalDownscale/5kASCII/Solar//2009/20090211.asc.gz"

#### This script will extract all 9 second cells from GRASS clear sky radiation that fall within a single .5 degree cell from AWAP realized radation files
#### The 9 second cells will be weighted by the mean values of all 400 cells within the .5 degree AWAP cell
#### These 9 second cell weighted values will then be multiplied by the .5 degree AWAP realized radiation
#### This will produce a 9 second resolution map of realized radiation

# Obtain command line arguments from .sh file

args=(commandArgs(TRUE))

# Evaluate the arguments for use in this script

for(i in 1:length(args)) 
	
	{
		
	eval(parse(text=args[[i]]))
 
	}
	
# Load libraries
library('SDMTools')

# Directories

dir.5k = '/home1/99/jc152199/MicroclimateStatisticalDownscale/AWAPASCII/RawAWAP/Solar/'
dir.250m = '/home1/99/jc152199/GRASS/ASCIIoutput/'
out.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/Realized_Radiation/'

# Read in an ASCII file size WTplusbuffer to use as an analysis mask

mask.250 = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/Analysis_Mask/Analysis_Mask_ASCII/AnalysisMask_WTplusbuffer_LatLong_WGS1984_250mres.asc')
#mask.250[,]=NA # Set all values to NA
mask.pos = as.data.frame(which(is.finite(mask.250), arr.ind=T)) #Get all row/col locations for the mask
mask.pos$lat = getXYcoords(mask.250)$x[mask.pos$row]
mask.pos$long = getXYcoords(mask.250)$y[mask.pos$col]

# Create a list of solar 5k files

solar.5k.files = list.files(dir.5k, pattern = '.grid.gz', recursive=TRUE, full.names=TRUE)
solar.5k.files = solar.5k.files[c(1:(length(solar.5k.files)-(365*5)-152))] # Select only files from 20040101 - 20090601

# The following loop creates a list of DOY's (i.e. 001-365) that is the same length as the list of solar.5k.files

DOYlist = NULL

for (yearx in c(unique(substr(solar.5k.files,97,100))))

	{
	
	t.DOYlist = solar.5k.files[(which(substr(solar.5k.files,97,100)==yearx))] # Grabs only files from a single year
	
	DOYS = sprintf('%03i',c(1:length(t.DOYlist))) # Creates a vector of numbers the same length as the number of solar files for that year
	
	DOYlist = c(DOYlist,DOYS) # Binds these vectors to a list and when finished this loop will create a vector of DOY's the same length as solar.5k.files
	
	}
	
	# Create a temporary copy of mask.pos
	
	t.pos=mask.pos
	
	# Select the appropriate DOY based on the position of the solar.5k.file
	
	t.DOY=DOYlist[which(solar.5k.files[]==tfile)]
	
	# Read in the Solar ASCII twice, one will have values converted to unique indices
	
	t.solar = read.asc.gz(paste(tfile,sep=''))
	t.solar.indices = read.asc.gz(paste(tfile,sep=''))
	
	# Assign each .5 degree cell from the AWAP solar file a unique index
	 
	tt = which(is.finite(t.solar.indices)) 
	t.solar.indices[tt] = c(1:length(tt)) 
	
	# Read in the appropriate 250 m Clear Sky Radiation File

	t.clearsky = read.asc(paste(dir.250m,'totalrad_',t.DOY,'.asc',sep=''))
	
	# Intersect mask positions with t.solar.indices to extract the unique index of that cell
	
	t.pos$index5k = extract.data(cbind(t.pos$lat,t.pos$long),t.solar.indices)
	
	# Intersect mask positions with t.solar to extract values for the 5k solar cells
	
	t.pos$value.5k = extract.data(cbind(t.pos$lat,t.pos$long),t.solar)
	
	# Intersect mask positions with t.clearsky to extract clear sky radiation values
	
	t.pos$clearsky = extract.data(cbind(t.pos$lat, t.pos$long), t.clearsky)
	
	# Aggregate values of clear sky radiation by 5k index value
	
	t.ag = aggregate(t.pos$clearsky, by=list(id=t.pos$index5k),FUN=mean)
	
	# Assign matching names to t.ag
	
	names(t.ag)[1]='index5k'
	names(t.ag)[2]='meanclearsky'
	
	# Create a column in t.pos to bind meanclearsky values into, populate it with NA's
	
	t.pos$meanclearsky=NA
	
	# Begin a for loop to match indices in the aggregate dataframe to the full t.pos dataframe, and populate the meanclearsky column with real values
	
	for (i in c(1:nrow(t.ag)))
	
		{
		
		t.pos$meanclearsky[which(t.pos$index5k==t.ag$index5k[i])]=t.ag$meanclearsky[i]
		
		cat('\n',i)
	
		}
	
	# Close loop	
	
	# Feed meanclearskyrad values back into mask.asc and write out with an appropriate label
	
	mask.250[cbind(t.pos$row,t.pos$col)]=t.pos[,8]
	
	write.asc.gz(mask.250,file=paste(out.dir,'realrad_',substr(tfile,97,104),'_WTplusbuffer_LatLong_WGS1984_9secres',sep=''))
	
	# Done
